function [A, B, Phi] = AIA_step1(I,Delta)
    [I_row, I_col] = size(I{1,1});
    A = zeros(I_row, I_col);B = A;Phi = A;
    for r = 1:I_row
        for c = 1:I_col
            S = AIA_step1_S(Delta);
            R = AIA_step1_R(I,Delta, r, c);
            X = S \ R;
            A(r,c) = X(1);
            Phi(r,c) = atan2(-X(3), X(2));
            B(r,c) = sqrt(X(2)^2+X(3)^2);
        end
    end
end